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Abstract. We investigate the influence of indirect connections, interre- 
gional distance and collective effects on the large-scale functional net- 
works of the human cortex. We study topologies of empirically derived 
resting state networks (RSNs) , extracted from fMRI data, and model dy- 
namics on the obtained networks. The RSNs are calculated from mean 
time-series of blood-oxygen-level-dependent (BOLD) activity of distinct 
cortical regions via Pearson correlation coefficients. We compare func- 
tional-connectivity networks of simulated BOLD activity as a function 
of coupling strength and correlation threshold. Neural network dynam- 
ics underpinning the BOLD signal fluctuations are modelled as excitable 
FitzHugh-Nagumo oscillators subject to uncorrelated white Gaussian 
noise and time-delayed interactions to account for the finite speed of 
the signal propagation along the axons. We discuss the functional con- 
nectivity of simulated BOLD activity in dependence on the signal speed 
and correlation threshold and compare it to the empirical data. 
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1 Introduction 

Despite important progress over past few years, little is known about brain 
functional connectivity (FC) at rest, i.e. under no stimulation and in the ab- 
sence of any overt-directed behaviour. In the studies on goal-directed mental 
activity, spontaneous brain activity has been considered as random enough to 
be averaged out across many trials [1] . However, well organized spatio-temporal 
low-frequency fluctuations (< 0.1 Hz) have been observed in blood-oxygen-level- 
dependent (BOLD) fMRI signals of a mammalian brain in the absence of any 
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stimulation or task related behaviour [2 3 4J . These well organized patterns of 
activity, suggest existence of underlying dynamics that governs intrinsic brain 
processes. 

Existing large-scale models of the intrinsic brain dynamics focus on the rela- 
tionship between functional and anatomical connectivity. They explore the range 
of conditions under which functional networks may emerge from known anatom- 
ical connections. In particular, roles of multiple time-scales in FC networks for- 
mation [516] . time-delays in the signal propagation between the network nodes 
and system noise |7l8j . local network oscillations [9110] and structural discon- 
nection [11] in the underlying dynamics have been investigated. Although RSNs 
reflect anatomical connectivity between brain areas comprising the networks in 
focus, FC cannot be understood in those terms alone [12113) . 

Here, we combine experimental and modelling approaches to investigate dy- 
namics underlying correlated behaviour of distant cortical regions. We choose to 
model the local node dynamics by excitable FitzHugh-Nagumo (FHN) neurons 
[14115] and use the resulting time-series to infer the neuronal activity at the 
network levels and subsequently low-frequency oscillations in the BOLD data 
[16] . We focus on identifying how global coupling strength and different network 
topologies affect connectivity patterns in simulated BOLD signals. 

The rest of this paper is organized as follows: In Sec. [2] we provide a detailed 
description of the methods used to generate the networks. In addition we intro- 
duce the model equations. This sets the stage for numerical simulations that are 
presented in Sec. [3j Finally we conclude with a discussion in Sec. [4] 

2 Methods and Models 

Subjects. We use resting state blood-oxygen-level-dependent (BOLD) fMRI sig- 
nals downloaded from the 1000 Functional Connectome Project website (http: 
|//www. nitric . org/[ ). 26 functional and anatomical scans from the Berlin study 
are considered in the analysis. Demographic data for the subjects participated in 
the study along with the technical details of the signal acquisition can be found 
at the Functional Connectome website. 

Data Preprocessing Steps in FSL. Image preprocessing is carried out 
using FSL (http://www.fmrib.ox.ac.uk) [17] . The preprocessing consists of 
the following steps: (a) temporal high-pass filtering (using Gaussian-weighted 
least-squares straight line fitting with sigma equals 100 s); (b) temporal low-pass 
filtering (using Gaussian filter with HWHM = 2.8 s); (c) slice-timing correction 
for interleaved acquisition (using Fourier-space time-series phase-shifting); (d) 
motion correction (using a six parameter affine transformation implemented in 
MCFLIRT); (e) spatial filtering (using Gaussian kernel of FWHM = 6 mm); 
(f) non-brain removal (BET brain extraction is applied to create a brain mask 
from the first volume in the fMRI data): (g) normalization to standard brain 
image (using 12 DOF linear affine transformation implemented in FLIRT each 
subject scan was transformed in MNI space - voxel size 2x2x2 mm). 
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Functional Connectivity. For the connectivity analysis we extract BOLD 
time-series from N = 64 cortical regions. These regions are adapted from stud- 
ies of functional segmentation of the human cortex |18ll9j . They consist of 30 
pairs of the inter-hemispheric homologues and 4 additional areas chosen along 
the midline. The full list of the cortical areas along with their anatomical (MNI) 
coordinates is provided in the work by Anderson et al. in Table 1 in the Supple- 
mentary material of Ref . [19] . We define the areas of interest as 5 mm spheres 
centered at the corresponding anatomical coordinates. For the given image res- 
olution of 2 x 2 x 2 mm voxel size (measured at 64 x 64 x 33 brain sites in 
total), there are 81 voxel time-series within each of the spherical areas. To ob- 
tain correlation matrix, describing connections between the regions, we calculate 
the linear correlation coefficient between any pair of the mean signal intensity 
of each of the 64 regions as: 

_ (V( Xl ,t)V(x 2 ,t)) - (V(x 1 ,t))(V(x 2 ,t)) 
r{XuX2) ~ a(V( Xl ,t))a(V(x 2 ,t)) ' [i) 

where V(x, t) represents activity of the region x at time t (averaged across 81 
voxels), a is the standard deviation, and (•) denotes temporal averages. The 
N x N correlations matrices are averaged across the ensemble of 26 subjects. 
The obtained correlation matrix {/y}, i,j — 1,...,N, is used to create FC 
maps between the cortical regions of interest. By definition, FC between two 
brain regions exists if their temporal correlation exceeds some predetermined 
value r, regardless of their anatomical connectivity. 

Thresholding: We obtain FC network topologies by thresholding the empir- 
ically derived FC matrix at different threshold values r. If fij > r, the corre- 
sponding element of the adjacency matrix is set to 1; otherwise it is set to 0. 
We simulate neural and BOLD activity only for r > 0.26, i.e. values for which 
connection density or topology cost in the brain networks is less than 0.5 [20] . 

Simulation of the Network Dynamics. To infer the BOLD signal in 
dependence on the network properties, we first simulate the underlying neural 
activity. We consider the neural dynamics on the network of N nodes or cortical 
regions, where local dynamics of the each node is represented by the homoge- 
neous FHN neurons. 

Simulation of the Neural Activity: Simulations of the neural network dy- 
namics take into account time-delays due to finite speed of signal propagation 
between the nodes as well as the presence of the system noise: 

N 

u t . = g[ui,Vi) - fijUjit - Atij) + n u (2) 

Vi = h(ui,Vi) + n v , (3) 

where c denotes global coupling parameter (c > 0), Ui and Vi are the activator 
and inhibitor variables of the i th neural population, /y are the elements of the 
FC matrix and Atij denote time-delays. n u and n v are two independent additive 
white Gaussian noise terms with zero mean, unity variance and noise strength 
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D. The functions g and h describe the local neural activity according to the 
FHN model, where we use the notation of Refs. |7|8j : 

u = g(u, v) = t + ju - y ^ (4) 
V = h(u, v) = (it — a + bv — I) , (5) 

T 

where / is magnitude of an external stimulus, which is assumed to be in our 
model of the resting state dynamics [7] . The system parameters are chosen as a = 
0.85, f3 = 0.2, 7 = 1.0 and r = 1.25 to render solutions with damped oscillatory 
behaviour of a node dynamics, i.e. at the onset of instability, in the absence of 
the connectivity in the network |7I8) . Representative time-series of the dynamics 
of an isolated node are shown in Fig. [l] for different noise strengths D = (no 
noise), D — 0.01 (small noise) and D = 0.05 (large noise). We solve the system 
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Fig. 1. Time evolution of the activator variable of the FitzHugh-Nagumo neural model 
for an isolated node under noise with different intensity. 



of coupled differential equations with time-delays and additive noise using the 
PYTHON-module PYDELAY [3T] . The algorithm is based on the Bogacki-Sampine 
method 22 23], which is also implemented in Matlab's DDE23. We calculate 
time-delays for a physiologically realistic value of the propagation velocity (v = 
7 m/s) via Atij = dij/v, where dij are approximated by the Euclidean distances 
between the centers of the spherical regions, i.e. network nodes i and j. A colour- 
coded representation of the d.^ values is shown in Fig. [2] We simulate 7.5 minutes 
of the real time using a time step of 1 ms. When time-delays, system noise and 
global coupling are taken into account, the simulated dynamics of the neural 
activity exhibit a behaviour that are exemplary shown in Fig. [3j 
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Simulation of the BOLD Activity: To infer the BOLD activity from the sim- 
ulated neural activity we assume the Baloon-Windkessel hemodynamic model 
[TB]- Simulated BOLD signals are band-passed in the frequency interval (0.01 - 
0.15) Hz to match low-frequency oscillations present in the experimental BOLD 
data and finally down-sampled to 2.3 s to match the temporal resolution of the 
MRI scanner. 



3 Results 

Empirical Functional Connectivity. FC matrices obtained from fMRI data 
and its binarized versions for different thresholds are shown in Fig. [4] (left) and 
the top panels Fig. [1] (right), respectively. In the latter figure, the threshold val- 
ues are chosen as r = 0.26, 0.38 and 0.5. r — 0.26 represents mean value of all 
correlations in the empiricaly derived FC and also the value at wich connection 
density (k) of the network drops below 0.5, (k — 0.48). Values r = 0.38 and 
0.5 are one and two standard deviations away from the mean, respectivelly. 




Fig. 2. Distance matrix {dij}, i,j 

Y between cortical regions in 

colour code, dij are calculated as Eu- 
clidean distance between centers of spher- 
ical regions in anatomical space. Matrix 
is ordered in a way that corresponding 
contra-lateral regions are symmetrically 
arranged with respect to the matrix cen- 
tre (right hemisphere regions, nodes 1 - 
30; left hemisphere regions, nodes 35 - 64). 
4 regions selected along the midline are 
placed in the middle of the matrix (nodes 
31 - 34). 
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Fig. 3. Time-series of neural network 
nodes dynamics for one highly and one 
sparsely connected nodes in visual cortex 
(VI, blue) and anterior intraparaetal sul- 
cus (AIS hIP2, red), respectively. Network 
dynamics are modelled with time-delays, 
large noise (D = 0.05) and weak cou- 
pling (c = 0.016). Parameters: a = 0.85, 
P = 0.2, 7 = 1.0 and r = 1.25. The inset 
shows the corresponding power spectra. 



Horizontal-plane view illustrate the anatomical maps of the corresponding net- 
works structure in the bottom panels of the Fig. Hi (right). In this representation 
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Fig. 4. (Left) Functional connectivity matrix constructed by calculating Pearson cor- 
relations on all pairwise combinations of the BOLD data from 64 cortical regions. The 
matrix is ordered as in the Fig. [2] The anti-diagonals reveal existing high correlations 
between contra-lateral regions. (Right) Top panels: Binarized functional-connectivity 
matrix at thresholds r = 0.26, 0.38, and 0.5. Each element is either black (if there is 
no significant connections between the regions) or white (if there is). Bottom Panel: 
Visualization of thresholded matrices in anatomical space by locating each region ac- 
cording to its x and y coordinates and drawing a link between significantly connected 
regions. 



each cortical region is treated as a network node with the connections or links to 
the other nodes in network if the matching correlation exceeds given threshold. 
FC networks show a pronounced left to right symmetry for different values of 
the correlation thresholds. However, distributions of the areas with most connec- 
tions depend on the threshold applied. The areas showing high numbers of links 
for the lower threshold network differ from the areas that appear as highly con- 
nected at higher thresholds. We expect that nodes with increased connectivity 
play an important role in the system's dynamics. 

Simulated Functional Connectivity. FC of the simulated BOLD activ- 
ity are shown in Fig. [5] The data are obtained varying correlation threshold r, 
i.e. network topology and global coupling strength c while keeping other model 
parameters at constant values. Correlations in the brain FC networks may not 
necessarily originate from direct connections. Therefore, we vary the network 
topologies by changing the correlation threshold r in the binary filter applied 
to the empirically derived FC network. We use the corresponding FC matrices 
to study the dynamics that result in the appearance of these correlations. Fur- 
thermore we uniformly scale all connections strengths c between the network 
nodes. For low correlation thresholds and weak coupling strengths (c < 0.05) 
the simulated BOLD signals are weakly correlated as the underlying neuronal 
activity does not show correlated behaviour either (data not shown). Increasing 
the coupling strength, positive correlations among the BOLD signals emerge in 
the simulated data and corresponding FC networks exhibit patterns of correlated 
activity similar to the empirically derived FC. 
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Fig. 5. Functional connectivity (FC) of simulated BOLD activity for different correla- 
tion thresholds r applied to the empirically derived FC matrix and different coupling 
strengths c. Other parameters as in Fig. [3] 



4 Discussion/Conclusion 



In this study we presented data obtained by combining both experimental and 
modelling approaches to explore dynamics underlying correlated behaviour of 
distant cortical regions. We showed that increasing global coupling strength be- 
tween the network nodes introduces correlations into simulated BOLD activity. 
The spatial patterns of the correlated activity resembled those observed in cor- 
responding experimental data. However, to get tangible measures of similarity 
between experimentally and theoretically derived FC networks our next step 
would be to apply graph-theoretical analysis to the obtained data. 

Our results are consistent with recent findings showing that stronger struc- 
tural couplings generates more globally connected and globally integrated BOLD 
signals . Since functional connectivity implies possible role of indirect anatom- 
ical connections further studies, showing how our findings relate to known anatom- 
ical connections between cortical regions could lead to the better understanding 
of FC networks formation. 
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